Existence of Positive Steady States for Mass 
Conserving and Mass- Action Chemical Reaction 
Networks with a Single Terminal- Linkage Class * 

b 

^ ' Santiago Akle^ Onkar Dalal^ Ronan M.T. Fleming-f 

^ , Michael Saunders^ Nicole Taheri'^ Yinyu Ye^ 

o 

o 

<N . 

Abstract 

We establish that mass conserving, single terminal-linkage networks 
of chemical reactions admit positive steady states regardless of network 
deficiency and the choice of reaction rate constants. This result holds 
O ' for closed systems without material exchange across the boundary, as 

. well as for open systems with material exchange at rates that satisfy a 

simple sufficient and necessary condition. 

Our proof uses a fixed point of a novel convex optimization formu- 
lation to find the steady state behavior of chemical reaction networks 
, that satisfy the law of mass-action kinetics. A fixed point iteration can 

be used to compute these steady states, and we show that it converges 
for weakly reversible homogeneous systems. We report the results of 
' our algorithm on numerical experiments. 

^ ; 1 Introduction 

One of the interests of systems biology is to deterministically model chemi- 
cal reaction networks. Models of such systems based on mass-action kinet- 
iJ> ■ ics depend on the kinetic parameters of the system. However, measuring 

^ ! these parameters experimentahy is difficult and error-prone. Thus, we seek 



> 



properties of chemical reaction networks that are independent of kinetic 
parameters. 
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In this work, we address the issue of existence of positive steady states, 
i.e., positive concentrations of species that will stay constant under the sys- 
tem's dynamics. In particular, we tackle the case in which a directed graph 
of chemical reactions forms a strongly connected component and the reac- 
tions conserve mass. We prove the existence of positive equilibria for such 
networks when they are closed systems, and extend our methods to open sys- 
tems where complexes are exchanged across the boundary at certain rates. 
How best to compute such a steady state remains uncertain; however, for 
closed systems we suggest a fixed point algorithm and provide results of our 
numerical experiments for several cases. 

1.1 Background 

Chemical Reaction Network Theory (CRNT), a mathematical theory for this 
problem, has its roots in the seminal work by Fritz Horn, Roy Jackson, and 
Martin Feinberg in O El El El [4] . We build on the notation used in these 
works and summarized by Gunawardena et al. in [6]. The system consists 
of a collection of species reacting collectively in some combination to give 
another combination of species in a network of chemical reactions. Let S 
be a set of m species and C be a set of n complexes. The relation between 
species and complexes can be written as a non- negative matrix Y € IR"*^"', 
where column yj represents complex j, and Yij is the multiplicity of species 
i in complex j. For example, the multiplicity of the species NaCl in the 
complex (H2O + 2NaCl) is 2. 

A reaction network is represented by the underlying weighted directed 
graph G{V,E), where each node in V represents a complex, each directed 
edge i ^ j denotes a reaction using i to generate j, and the positive edge 
weight ki-^j is the reaction rate. The matrix A G IR*^^" is the weighted 
adjacency matrix of the graph, where Aij = ki^j. Define D := diag(74l), 
where 1 is the vector of all ones, and := A'^ — D. The elements of this 
new matrix will be {Ak)ij = kj^i for i / j, and {Ak)ii = — kj^i, so that 
Ajl = 0. 

Let c G be the vector of concentrations of each species, and b € 
IR"^ the vector of species exchange rates across the network boundary. We 
establish necessary and sufficient conditions on the external exchange rates 
b in order to determine the existence of steady state concentrations, and to 
be able to compute it given a specific b. 

We define Tp{c) : IR™ — )• IR" to be a nonlinear function that captures 
mass- action kinetics: 

V'.(c)=ncr"- 

i 

The change in concentration over time can be described by the system of 
ordinary differential equations 

c = YAkiPic) - b. 
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Hence, steady state concentrations for a chemical reaction network are any 
non-negative vector c* € R'" such that YAkiJj{c*) = b. Equivalently, a pair 
{c*,v*) with c* G R"* and v* G wih be a steady state if it satisfies the 
conditions 

YAkV* = b (FB) 
i;{c*) = V* (MA) 

or flux-balance and mass-action, respectively. Thus, finding steady state 
concentrations is equivalent to finding a vector v* that satisfies both (IFBh 
and pyfA]) for some vector c* € R™'. 

Observe that if (c*, v*) is a positive steady state, then from the definition 

of ip{c*), 

log(c*) = log(V'(c'')) = log(?;''). (MA-log) 

We refer to this alternative condition as the logarithmic form of the mass- 
action condition ()MAp . In this notation, systems with the property of mass 
conservation can be characterized by the following definition. 

Definition 1. A chemical reaction network is mass conserving if and only 
if there exists a positive vector e G R™ such that 

e^YAk = 0, (1) 

where e denotes the molecular weights of the species (or atomic weights if 
the species are elements). 

The connectedness of the networks is captured in the following definition 
of a terminal-linkage class. 

Definition 2. A terminal-linkage class is defined as a set of complexes C 
such that for any pair of complexes {i,j) G C there exists a directed path in 
the graph G that leads from i to j. 

We further restrict our analysis to a class of weakly reversible networks. 

Definition 3. A chemical reaction network is weakly reversible if it is 
formed exclusively by one or more terminal-linkage classes. 

A reaction network that consists of exactly one terminal-linkage class is 
called a single terminal-linkage network. Reversibility, at least in a weak 
sense, is a prerequisite for steady states with positive concentrations for all 
species, as suggested by simple examples like a single non-reversible reaction. 

Next we define a stoichiometric subspace and deficiency for a network. 

Definition 4. A stoichiometric subspace S is the subspace defined by the 
span of vectors yj — yi, where yj,yi are the columns ofY representing com- 
plexes i and j, for each reaction pair i ^ j in the network. 
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Definition 5. The deficiency of a network is defined as 5 = n — t — s, 
where t is the number of terminal-linkage classes and s is the dimension of 
the stoichiometric subspace, also known as the stoichiometric compatibility 
class. 

In one of tlieir early works, Horn and Jackson |9j analyzed mass-action 
kinetics for closed systems (with 6 = 0) and defined a class of equilibrium 
points called complex-balanced equilibria, and defined systems admitting 
such an equilibrium to be complex-balanced systems. These closed systems 
are shown to satisfy the quasi-thermo static and quasi-thermodynamic con- 
ditions regardless of the kinetic rate constants. Following this, Horn [7j 
also proved necessary and sufficient conditions for existence of a complex- 
balanced equilibrium. In [8], Feinberg and Horn used the existence of a 
Lyapunov function to show the uniqueness of the positive steady state in 
each stoichiometric compatibility class, which is equivalent to specifying all 
the conserved quantities of a system. Later Feinberg [3l proved two the- 
orems, now famously known as Deficiency 0-1 theorems, that provide the 
analysis of positive steady states for a class of networks with deficiency 
or networks with deficiency 1 but with each terminal-linkage class having 
deficiency less than 1. For this restricted class of closed systems, the ex- 
istence of a positive steady state is given by Perron-Frobenius theory for 
a positive eigenvector. Other work in this area is from the perspective of 
dynamical systems and aimed toward proving two open conjectures: Global 
Attractor Conjecture and Persistence Conjecture [1]. Another approach us- 
ing parametrized convex optimization to compute a non-equilibrium steady 
state is given in [5]. 

To the best of our knowledge, the vast majority of CRNT research stud- 
ies closed complex-balanced systems, which, by definition, admit a complex- 
balanced equilibrium. In the notation above, a complex-balanced equilibrium 
exists when the vector of concentrations c satisfies Akip{c) = 0, i.e., the vec- 
tor V'(c) belongs to the null space of A^. It can be shown that a network with 
linearly independent complexes will have deficiency 5 = 0. However, if some 
of the complexes are linearly dependent (as shown in the example in Section 
3), there are systems that are not complex-balanced yet admit concentra- 
tions in equilibria where Ai:ip{c) is in the null space of Y and Akipic) 7^ 0. 
Though the condition of complex-balance is sufficient for thermodynamic 
consistency, p] shows that it is not necessary. Also, for open systems with 
material exchange across the boundary, complex-balance is not defined. In 
order to handle open systems, these works hint at extending the system us- 
ing a pseudo 0-complex and adding pseudo reactions. However, it is unclear 
how to choose the pseudo kinetic rates such that the positive eigenvector so- 
lution of the extended system will achieve the given external exchange rates 
b. From the point of view of systems biology and bio-chemical engineering, 
analyzing the behavior of a cell under different exchange conditions b is very 
important to control and engineer the cell, for example studying the desired 
effects in pharmacology, or producing specific metabolites in bioreactors. 
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In this paper, we extend the previous work on two accounts: 1) we prove 
existence of positive equihbria in closed systems for some reaction networks 
that do not satisfy the necessary conditions of the Deficiency 0-1 theorems 
(are not necessarily complex-balanced), and 2) we provide a necessary and 
sufficient condition on the external exchange rate b for some open systems to 
admit a positive steady state. We use a fixed point of a convex optimization 
problem, with an objective function similar to the Helmholtz function de- 
fined in The fixed point of this mapping gives the required steady state. 
We prove the existence of a positive steady state for any weakly reversible 
chemical reaction network with a single terminal-linkage class. We strongly 
believe that this can be extended to systems with multiple terminal-linkage 
classes, as supported by our computational results for randomly generated 
networks. Section [32] gives a detailed analysis of a toy network to emphasize 
this claim. 



2 A Fixed Point Model 

Our main result establishes that for any set of positive reaction rates k G R" 
and any b in the range of YAk , a single terminal-linkage network will admit a 
positive solution pair {c,v) that satisfies the laws (IFBp and (]MA|1 . We show 
this by defining a positive fixed point of a convex optimization problem, and 
establishing an equivalence between the positive fixed point and positive 
solution to the equations. 

We construct a fixed point mapping of a linearly constrained optimiza- 
tion problem such that the logarithmic form of the mass-action equation 



(MA-log) is an optimality condition, and hence any solution to this opti- 



mization problem will also satisfy (MA-log) 



To define the mapping, let b = YAkrj for some rj € R" and observe that 
for arbitrary s € R" we can write b = YA'^{ri + s) — YD{r] -\- D~^A^s). In 
particular, we choose s positive and large enough so that r]^ := r] -\- s and 
r]~ := r/ -|- D~^A^s are both positive. Also, from Definition [H e'^b = and 
thus 

e^YDr]- = e'^YA^ri+. (2) 

Define /u := (r, tq) G IR''"^"^ to be a vector parameter. Observe that if the 
parametric convex optimization problem 

minimize v"^ D{log{v) — 1) + fo(logfo — 1) 

subject to YDv + YA^rj+vo = YA^r + YDrj^r^ : y (3) 
{v,vo) > 

has a positive solution {v*{fj,),VQ{fj,)), then the optimality conditions 
YDv*{n) + YA^r]+v^ = YA^r + YD-q-r^ 

(yAV)V(^) =logK(//)) 
{v*{fi),v*o{fi)) >0 
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are well defined. Since D is nonsingular, the second optimality condition 
is equivalent to (MA-log), for c*(/x) := e^**^^^, where the exponent is taken 
element- wise. Hence, the equation ( |MA-log ) holds and c*{fi) satisfies mass- 
action. We show that for some parameter fi, the equality 

YAkV* = v*ob (5) 

is satisfied. This implies that for 6 = 0, both (IFBp and (MA-log) are satisfied, 
and the solution is attained. For the case where 6 7^ 0, we can construct a 
corresponding solution so that (IFBh holds. 

Note that the nonlinear program ([3]) is strictly convex, so for any feasible 
(r, ro) there is a unique minimizer. That is, the mapping 

(r,ro)^(t;^(/x),t;S(/i)) (6) 

is well defined. If /i = (r, ro) is a fixed point of ([6]), then the linear equality 
constraint in ([3]) implies 

YDv*{fi) + YA^rj+v^ifi) = YA^v*{fi) + YDr]-v^{fi) 

or, equivalent ly, 

YAkV^fi) = - D)v*{fi) = v*{YA'^rj+ - YDrj-) = v*{il)b. 

Therefore, if such a fixed point exists, the solution v*{jl) at this fixed point 
will satisfy ([5]). For simplicity, we henceforth refer to the optimal solution 
variables {v*{fj,),VQ{fi)),y*{fi) as {v*,VQ),y*, but acknowledge their depen- 
dence on fi. 

Theorem 1. For any mass conserving, mass-action chemical reaction net- 
work and any choice of rate constants k > 0, there exist nontrivial fixed 
points for the mapping ([6]) . 

Proof. Brouwer's fixed point theorem states that any continuous mapping 
from a convex and compact subset of a Euclidean space O to itself must 
have at least one fixed point. 

Let {v*,Vq) be defined as in ([6]) and let 7 be a positive fixed scalar. 
Define the set 

n = {{v, vq) e 5R"+^ : (u, vq) > 0, e^YDv + e'^YA^r]+vo = 7} , 

where e is defined in ([1]). According to Brouwer's fixed point theorem, if 
the parameter (r, rg) G O ensures that the corresponding solution to the 
optimization problem {v*,Vq) G Q, then there is a fixed point such that 
the parameter and the solution are equal, i.e., there exists a /x such that 
= (r,ro) = {v*{fi),v^{fi)). 
The set is bounded and formed by an intersection of closed convex 
sets, and hence is convex and compact. Moreover, the mapping ^ (v* ,Vq) 
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is continuous. Since problem ^ is feasible for any /i € il, the mapping 
Q 3 fj, {v*,Vq) is well defined. 

To show that the image of O under the mapping (r, ro) — >■ (v*,Vq) is in 
0, first observe that by the bounds in (l3|, {v*,Vq) > 0. Using the equality 
constraints, Definition ([T]) and Equation ([2]), we have 

e^YDv* + e^YA^r]+v^ = e^YA^r + e^YDr]~ro 

= e^YDr + e^YA^rj'^ro = 7, 

and thus {v*,Vq) G 17. 

Therefore, under the mapping (r, tq) {v*,Vq), (r, tq) € implies 
(f*,f;Q) G 0, and the mapping must have a fixed point. Moreover, since Q 
does not contain the zero vector, the fixed point (s) are nontrivial. 

□ 

Note that the value of YDv* + YA'^'q^VQ is the rate of consumption of 
each chemical species and YA^v* + YDr]~VQ is the rate of production of 
each chemical species. At the fixed point, the equality YDv* + Y A^t^'^Vq = 
YA'^v* + YDri~VQ defines a steady state. The set Q. defines the parameter 
7 = e^(YDv* + YA^tj'^Vq); since the vector e can be interpreted as an 
assignment of relative mass to the species, 7 can be interpreted as the total 
amount of mass that reacts per unit time at the steady state. Therefore, 
looking for fixed points in il. corresponds to looking for steady states where 
the amount of mass that reacts in the system is prescribed. 

We have established the existence of a nontrivial fixed point fi of the 
mapping 9 // ^ (v*,Vq) € fi. Moreover, we have shown that when 
the associated minimizer {v*,Vq) is positive, it is a solution to ()MAp and 
to YAkV* = VqB. However, in the case when some entries of v* are zero, 
the objective function of ([3]) is non-differentiable and we cannot use the 
optimality conditions to show that ()MAp holds. 

2.1 Positive fixed points in single terminal-linkage networks 

We now consider the case when the network is formed by a single terminal- 
linkage class and show that if // is a fixed point of the mapping ([6]), the 
minimizer (//)), and therefore //, is positive. 

Lemma [U shows that if problem (j3]) has a feasible point with support 
J, the minimizer {v*,Vq) will have support at least J. Lemma [2] uses the 
single terminal-linkage class hypothesis to show that at a fixed point, there 
is a positive feasible point. These two Lemmas imply that at a fixed point 
fi, the minimizer will be positive. Finally Theorem [2] shows that if -Oo 7^ 1 
at the solution, we can construct another solution for which vq = 1. This 
establishes that there is a nontrivial steady state for the network. 

To complete the argument we must prove Lemmas [H [2] and Theorem [2j 

Lemma 1. The support of any feasible point of Problem ^ is a subset of 
the support of the minimizer {v*,Vq). 
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Proof. Let v G 5}?"^ be any of the feasible points with the largest support 
and let z be any feasible direction at v. By construction, for all a in some 
interval [£, u] the points Va := v + az are non-negative and feasible. The 
interval can be chosen so that when a = I and when a = u, one new bound 
constraint becomes active. This implies that supp(t;£) and supp(fu) are 
strictly contained in supp('5), and supp(?jQ,) = supp('u) for a G {i,u). 

Without loss of generality, we assume i < < u, since i and u will not 
be of the same sign; if ^ = and u > 0, any point Va can be written as a 
convex combination of v and v + uz, and thus has support as large as v. 

Define the univariate function 

g{a) := (p{v + az), (7) 

where cp is the objective function of ([3]). We will establish that as a — > Z 
the derivative g'{a) — oo, and as a ^ n the derivative g'{a) — >■ oo. Thus, 
by the mean value theorem, there must exist a zero of the function g in 
the interior of the interval [l,u]. Since this function is strictly convex, if a 
stationary point exists in the interior of the interval, the function value at 
the stationary point must be smaller than at the boundary. 

Observe that if we let dj, for i € [1, . . . n], be the diagonal entries of D 
and dn+i = 1, we can write 

n+l 

g{a) = '^{v + OiZi)di \og{vi + azi). 

i=l 

An important observation is that if some entry Vj = then Zj = 0, 
otherwise Va would have a larger support for some a ^ 0. This implies that 
{va)j = for all entries where Vj = 0. If we let J be the set of nonzero 
entries of and L be the subset of J formed by the entries that tend to 
zero as a ^ then 

^'(a) = X] Zidi{\og{vi + azi)) 

= ^ Zidi{\og{vi + azi)) + ^^Zidi{\og{vi + azi)). 

As a ^ the first summation will approach a finite value. Since > for 
all i G -L, the entries in the logarithm of the second sum tend to zero and 
the term will diverge to — oo. 

Similarly, let U be the subset of J formed by the entries that tend to 
zero as a ^ u. Observe that for these entries, < and 

g'{a) = ^ Zidi{\og {vi + azi)) + ^ Zidi{\og [vi + azi)). 
The first sum will tend to a finite value and the second will diverge to oo. 
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Now, assume that for some there is a feasible point {vjVq) with larger 
support than the minimizer (v* ,Vq) of problem ([3|). Since (f *, Vq) has smaller 
support, we can write {v*,Vq) = {v,vo) + a*z where a* is on the boundary 
of the corresponding feasible interval. By the previous argument, there is a 
value of a 7^ a* in the interior of the interval such that {v*, Vq) = {v, VQ) + az 
has a lower function value than {v*,Vq), contradicting its optimality. 

Therefore, by the mean value theorm, there must exist a stationary point 
of g strictly in the interior of the interval [l,u] at which the function value 
is smaller than at the boundary. Moreover, the optimal point will have at 
least the support of any feasible point. 

□ 

Lemma 2. If the network is formed by a single terminal-linkage class, when 
Problem ^ is parametrized by a fixed point fi, there exists a positive feasible 
point {v,vo). 

Proof. Let Problem ([3]) be parametrized with a fixed point fi, and let {v,vo) 
be both the minimizer and the fixed point. We prove by contradiction that 
no entry of the minimizer {v,vo) can be zero. Observe that by the definition 
of 17 the origin is not contained in the set, and therefore the fixed point 
cannot be identically zero. 

First, assume that vq > and observe that p := {D^^A^v + r]~VQ,0) is 
a feasible point. Since r]~ was chosen to be positive and D~^A^ has no zero 
columns, the support of p are the first n entries of the vector. A convex 
combination of (£','0o) and p will be feasible and have full support. 

Now, assume that = and some entry of v is nonzero, and observe that 
{D-^A^v,Qi) is feasible. A convex combination of {D ^A^v,{)) and (?),0) is 
feasible and its support contains the union of the supports of the two vectors. 
That is, for (3 E [0,1] the point {v,vo) = p{D-^A^v,Q) + (1 - /3)(t},0) is 
feasible, and using the fact that the support of D~^A^v is the support of 
A^v along with Lemma [H 

(supp(i), 0) U supp(^^i), 0)) C supp(i), 0). (8) 

This relation can be used inductively to show that there is a feasible point 
with support at least as large as the union of the supports of {{A^Yv, 0) for 
all positive powers of p. 

The single terminal-linkage class hypothesis implies that for any pair 
(i, j) G [1, . . . , n] X [1, . . . , n], there exists a power p large enough such that 
{A^Y-- > 0. More importantly, if Vj > 0, then there exists a p such that 
[{A^yvj]i > for all i. Therefore, if = and v ^ 0, there is a feasible 
point {v, 0) such that v > 0. 

Finally, if v > then there is a scalar < a small enough such that 

0<v- D-^A^7]+a, 



9 



and then the equahty 

YD{v - D-^A^T]+a) + YA^r]+a = YDv = YA^v 

imphes that the positive point {y — D^^A^rj^a, a) is feasible. 

Therefore, if a network is formed by a single terminal-linkage class, the 
problem ([3]) has a positive feasible point {v^vq). □ 

Theorem 2. For a mass conserving single terminal-linkage network, there 
exists a concentration c > such that Y Akip[c) = b if and only ifb is in the 
range of YAj. ■ 

Proof. We have shown that there exist positive vectors c € R™', v € R" 
such that Y'^ log(c) = log(v) and YAf^v = vob. In other words, we have 
proven that there is a positive vector c and a positive scalar a such that 
YAkip{c) = ab. If we can construct a new concentration vector c > that 
satisfies V'(c) = ^ipic), then 

YAk^P{c) = -YAki^ic) = b 
a 

and the steady state concentration c satisfies (|FBp and (|MAp . 

First, we argue that the vector of all ones, 1 G IR,"', is in the range of 
when the network consists of a single terminal-linkage class and is mass 
conserving. The condition of mass conservation implies that e'^YAf^ = 0, 
or equivalently, Y'^e G 7\A(j4^). Since A^ is the Laplacian matrix of a 
strongly connected graph, M^A/^) = {f31 : j3 £ R}; thus, for some value /3, 
Y'^e = (31. 

Now, observe that log (^V'(c)) = log(V'(c)) — llog(a) = y-^log(c) — 
llog(a), where log(c) is an entry-wise logarithm of the vector and log(a) 
the scalar logarithm. Moreover, since 1 is in the range of Y"^ , say Y'^6 = 1 
for some 6 € R'", then log(a)l = Y^{log{a)6). Thus, if we define c as the 
vector that satisfies log(c) = log(c) — log{a)6 , then 

y^log(c) = y^log(c) - log(a)l, 

which implies that 

V'(c) = -ipic). 
a 

Therefore, the inhomogeneous system has a solution if the underlying graph 
of the network is formed by a single terminal-linkage class and the network 
is mass conserving, regardless of the kinetic parameters. 

□ 

3 Numerical Experiments 

The results in this section depend on the calculation of positive steady states 
for weakly reversible networks. We use the algorithm described below to 
solve for the associated fixed points. 
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3.1 Numerical method for finding fixed points 



Given an initial positive point v G R" and a small tolerance r, we use the 
following fixed point iteration, Algorithm [H to find a parameter fj, = (r, tq) 
to the problem ([3]) such that {v*,Vq) = (r, tq). 

Algorithm 1 Fixed Point Iteration to Find a Steady-State Concentration 

1: {v\v*q) ^ {V,0) 

2: (r,ro) ^ 

3: while \\YAkV* - (YAkr]~^ + YDr]^)vQ\\oo > r do 
4: {v*,Vq) ^ unique solution of ([3]) 
5: (r,ro) ^ i(r,ro) + i(^;^^;S) 

6: end w^hile. 



Step [4] in the while loop of Algorithm [T] requires solving the linearly 
constrained convex optimization problem Our implementation uses the 
PDCO package [13] to solve this problem. 

Provided that at each iteration k, the unique solution of ^ satisfies 
v*{nk) > and the minimization is solved with sufficient accuracy, the 
optimality conditions for ^ will imply that for all /c = 1, 2, 3, . . ., 



for some small value of e, where A* is the Lagrange multiplier of the linear 
equality constraint at the solution that corresponds to the logarithm of the 
concentrations. Thus, if the iteration converges to a fixed point (r, rg) = 
{v*,Vq), then at this fixed point (jMAj) will be satisfied to a precision e and 
()FBp will be satisfied to a precision r. 

Algorithm [1] has been tested extensively on randomly generated net- 
works with noteworthy success. The results of our experiments are shown 
in Section 13.41 

3.2 An example network 



Figure 1: Example network with two terminal-linkage classes 

In this section we consider the toy network shown in Figure [TJ The 
number of complexes is n = 7, the number of terminal- linkage classes is / = 2, 
and the stoichiometric subspace S = span{j4 + E — C,C — A — D, B — C} has 



y^A*-log(^;*(^fc))||oo <e 
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dimension s = 3. Therefore, the deficiency of this network is (5 = 7—2—3 = 2, 
and hence neither of the Deficiency 0-1 theorems [3l H] can be applied to 
calculate equilibrium points. 

However, since this network is weakly reversible, intuition suggest that a 
non-zero steady state exists. We use Algorithm [T] to solve for the fixed point 
described in Section [21 obtaining a positive steady state. Figure [2] illustrates 
the convergence of the fixed point iterations to the steady state. 




Figure 2: Concentration convergence with iteration 

Figure [3] illustrates the change in steady state as a function of the to- 
tal mass in the system, where total mass is defined as 7 = e^(YDv* + 
YA^tj'^Vq), as described in Theorem [TJ The experiment shows that as the 
total mass increases, species A, B and C adjust linearly to the additional 
mass, while species D and E stay at the same levels. This linear growth in 
species A, B and C can be explained analytically by the fact that the vector 
1 lies in range of Y'^ . 

3.3 Generating suitable networks 

This section describes the sampling scheme used to generate random mass 
conserving chemical reaction networks with a prescribed number of terminal- 
linkage classes. The output of the method is a network with n complexes, m 
species, and £ strongly connected components, where i is the desired number 
of terminal- linkage classes. 
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Figure 3: Equilibrium dependence on total mass 

First, we iteratively generate Erdos-Reyni graph^ with m nodes until 
we sample a graph with £ strongly connected components; call this graph 
G(y , E). We give each edge in £^ a weight of an independent and uniformly 
distributed value in the range (0,10]. These edge weights represent the 
reaction rates between complexes. 

To generate the stoichiometry, we define a parameter r as the maximum 
number of species in each complex. Each complex is constructed with a 
random sample of q species, where 5 is a random integer in [1, r]. All samples 
are done uniformly and independently. Finally, we assign the multiplicity of 
each species in a complex with independent samples of the absolute value of 
a standard normal unit variance distribution. To ensure mass is conserved, 
we normalize the sum of the stoichiometry of the species that participate in 
a complex to one, so that = 1 and = 0. 

3.4 Convergence of the Fixed Point Algorithm 

This section illustrates the convergence of Algorithm [1] on large networks 
that consist of either a single terminal-linkage class or multiple terminal- 
linkage classes, i.e., weakly reversible networks. 

Algorithm [T] produces sequences that, up to a small tolerance, satisfy 

^An Erdos-Reyni graph is a directed unweighted graph. Each edge is included with 
probabiUty p and all edges are sampled iid. 
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pvlA]) at every iteration. Ideally, the infeasibility with respect to ()FBp also 
monotonically decreases until convergence. Our extensive numerical exper- 
iments indicate that this is in fact the behavior for homogeneous networks 
that are weakly reversible. 

Mass balance infeasibility vs iteration count for single linkage network 
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Figure 4: Typical Infeasibility of (IFBp vs. Iteration, for network with a 
single terminal-linkage class 

Mass balance infeasibility vs iteration count for multiple linkage network 
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Figure 5: Typical Reduction in Infeasibility of ()FBp . for network with two 
terminal-linkage classes 

Figure [4] displays the sequence of the infeasibilities ||y^A:Wfc||oo at each 
iteration k in Algorithm [U for a network with a single terminal- linkage class, 
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50 species and 500 complexes, where at most 10 species participate in each 
complex. Figure [5] displays the analogous sequence for a network of equal 
size and two terminal-linkage classes. We have observed this (apparently 
linear) convergence rate consistently over all generated networks, regardless 
of the number of terminal-linkage classes. 
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Figure 6: Average number of iterations for single terminal-linkage class net- 
works 

We have also investigated the number of iterations necessary for Algo- 
rithm [1] to converge on networks of different sizes, with either one or two 
terminal-linkage classes. Figures [6] and [7] display the mean number of it- 
erations necessary for convergence on networks ranging from 100 to 5000 
complexes, where each average is taken over 20 instances per network size. 
Notably, the average number of iterations increases less than 10% as the 
network size grows fifty- fold. 

In future work, we plan to prove theoretical results on the existence of 
positive equilibria for chemical reaction networks with multiple terminal- 
linkage classes. However, our comprehensive numerical experiments seem to 
indicate that even for networks with more than one terminal-linkage class, 
there exists at least one positive fixed point of problem ([3]), and the iterates 
of Algorithm [1] converge to such a fixed point. 

Finally, while conducting this work, we were made aware of related work 
by Deng et al. [2]. Their work extends that of Feinberg and Horn by proving 
that weak reversibility is a necessary and sufficient condition for existence 
of positive equilibria. While their proof is for closed systems (with 6 = 0), 
it is not clear how to use the 0- complex in their extension to solve for any 
admissible b, as mentioned in Section 11.11 More importantly, our proof is 
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Figure 7: Average number of iterations for networks with two terminal- 
linkage classes 



based on a less complicated convex optimization formulation, which gives a 
method to calculate numerical solutions using a fixed point algorithm. 

Acknowledgment: We gratefully acknowledge Anne Shiu and her student 
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